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ABSTRACT 

The cosmic ray current-driven (CRCD) instability, predicted by|Bell (2004), consists of non-resonant, 
growing plasma waves driven by the electric current of cosmic rays (CRs) that stream along the 
magnetic field ahead of both relativist ic and non-relativistic shocks. Combining an analytic, kinetic 
model with one-, two-, and three-dimensional particle-in-cell simulations, we confirm the existence of 
this instability in the kinetic regime and determine its saturation mechanisms. In the linear regime, 
we show that, if the background plasma is well magnetized, the CRCD waves grow exponentially at 
the rates and wavelengths predicted by the analytic dispersion relation. The magnetization condition 
implies that the growth rate of the instability is much smaller than the ion cyclotron frequency. 
As the instability becomes non-linear, significant turbulence forms in the plasma. This turbulence 
reduces the growth rate of the field and damps the shortest wavelength modes, making the dominant 
wavelength, A^^, grow proportional to the square of the field. At constant CR current, we find that 
plasma acceleration along the motion of CRs saturates the instability at the magnetic field level such 
that va ^ Vd^cr^ where va is the Alfven velocity in the amplified field, and Vd^cr is the drift velocity 
of CRs. The instability can also saturate earlier if CRs get strongly defiected by the amplified field, 
which happens when their Larmor radii get close to A^^. We apply these results to the case of CRs 
propagating in the upstream medium of the forward shock in supernova remnants. If we consider only 
the most energetic CRs that escape from the shock, we obtain that the field amplification factor of 
~ 10 can be reached. This confirms the CRCD instability as a potentially important component of 
magnetic amplification process in astrophysical shock environments. 

Subject headings: ISM: magnetic filed — cosmic rays — supernova remnants — jets and outfiows 



1. INTRODUCTION 

High energy particle acceleration and magnetic field 
amplification appear to be tightly related phenomena in 
many astrophysical environments. For instance. X-ray 
observations of synchrotron emission from ultrarelativis- 
tic electrons in young sup ernova remnants (SNRs) (Volk^ 
et al.||2005[ |Ballet||2006l [Uchiyama et al.||2007[ ) suggest 
that electrons are etticiently accelerated in these environ- 
ments, and that the ambient magnetic field in the dow- 
stream medium of SNR forward shocks is amplified by a 
factor of ~ 100 compared to its typical value in the inter- 
stellar medium. Also, cosmic rays (CRs) with energies up 
to rsj lO^^eV are believed to originate in SNRs. Calcu- 
lations based on the di f fusive shock acceleration mech- 
anism p5ym sky"1977^ 'Axford e t al]|1977| |Bell||1978 



B landlord ^ O striker 1978) show that such high ener- 



gies can only be reached if CRs are efficie ntly confined 
to the remnant (Lagage & Cesarsky 1983), a condition 
that would be eased by a substantial magnetic field am- 
plification. 

While it is supported by direct observations of SNRs, 
the process by wh ich the field amplification takes place 
is still a mystery. 
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( 2004 ) suggested the possibility 
of magnetic ampliiication driven by non-resonant CRs 
propagating in the upstream region of shocks. This in- 
stability is different from the Alfven w ave amp lification 
due to cyclotron resonance wi th CRs (Kulsrud & Pearce 



1969[ [McKenzie fc yolk||1982[ ). Instead, it only requires 



positively charged CRs propagating along a background 
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magnetic field, Bq^ with Larmor radii much larger than 
the wavelength of the wave. This condition is expected 
to be satisfied by diffusively shock-accelerated CRs in the 
upstream region of both relativistic and non-relativistic 
shocks. While individual CRs are relativistic, on average 
they stream with respect to the upstream plasma with 
a drift velocity, Vd,cr^ that depends on the distance from 
the shock. The lowest energy CRs are efficiently con- 
fined to the shock vicinity by the upstream turbulence, 
and drift with the shock with Vd,cr ^ Vsh^ where Vsh is 
the shock speed. The higher energy CRs are less affected 
by the upstream scattering, so they tend to escape more 
easily from the shock. Thus, as the distance from the 
shock increases, Vd^cr will increase, asymptotically ap- 
proaching c/2 for non-relativistic shocks (for relativistic 
shocks, Vd^cr ^ c everywhere). Also, considering that 
electrons are more affected by radiative losses than ions, 
it is reasonable to think that at some distance from the 
shock CRs will be mainly positively charged particles. 
Their drift will drive a constant current, Jcr, through 
the upstream plasma. This current will be compensated 
by the opposite return current, J^et ^ —Jcr^ provided by 
the background plasma. If there is a small magnetic field 
perturbation Btr^ perpendicular to 5o, a Jret x ^tr force 
will push the background plasma transversely. A force 
of the same magnitude will push the CRs in the opposite 
direction, so that the net force acting on the background 
plasma-CRs system vanishes. However, given the high 
rigidity of the CRs, only the background plasma will ex- 
perience a significant transverse motion. For a helical 
magnetic perturbation Btr^ this transverse motion will 
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stretch the magnetic field fines, producing a n ampfifica- 
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tion of Bfr- Using analytic MHD analysis, 
sfiowed tfiat for a rigfit-fianded circularly po 
tromagnetic wave, tfiis amplification would be exponen- 
tial, and faster than the resonant instability^. 

The linear dispersion relation of this cosmic ray 
current-driven (C RCD) ins tability has been calculated 
using bot h MH D CBeir2004) and kinetic treatments ( Re- 
pille et al 



2006 Amato&Blasi2008). These works found 



that, if Jcr is kept constant, the instability will grow at 
a prefered wavelength Xmax = Bqc/Jcv^ with a growth 
rate ^rnax = Jcr{j^ 1 9(?Y^'^ ^ whcrc p is the mass density 
of the background plasma. 

The non-linear evolution of the CRCD instability has 
been studied making use of both MHD and particle-in- 
cell (PIC) simulations. The MHD studies |BeIT|[2QQ4l 
|2005t [Ziraka shvili et aL] 2008 ) have shown a substantial 
amplification of the ambient magnetic field and the for- 
mation of turbulence, which is characterized by promi- 
nent density fiuctuations in the plasma. They have estab- 
lished a saturation criterion that depends on the wave- 
length, A, of each mode and that is given by Bgati^) ^ 
B^X/Xmax- This saturation criterion would imply that, 
if Jcr is constant, the field never stops growing but only 
migrates into longer wavelengths. This migration would 
be such that the dominant wavelength, A^^, goes roughly 
as Xd ~ Xmax{B/Bo), where B is the mean value of the 
magnetic field (Bell 2004). The saturation in MHD sim- 
ulations with constant C R current is either not observed 

J Zirakashvili et al.|2008 ) or is due to the s ize of magnetic 
uctuations reaching the size of the box ([BeII|[2004). 
One important limitation of the MHD simulations is 
that they cannot follow the instability at arbitrarily low 
densities, which makes it difficult to model large plasma 
density fiuctuations properly. Also, the MHD simula- 
tions do not include the back-reaction on the CRs, which 
must play a fundamental role in the saturation of the in- 
stability. Both difficulties can be potentially resolved by 
fully kinetic PIC simulations. 

A first atte mpt to use PIC s imula tions for this problem 
was made by Niemiec et al.l (2008). Even though their 
results show niagnetic amplification, it accurs at a signif- 
icantly lower rate and through a kind of turbulence that 
essentially differs fr om the circ ularly polarized, growing 
waves predicted by Bell (2004). This raised a question 
about the existence of the CRCD instability beyond the 
MHD approximation. 

In this work we confirm the existence of the CRCD 
instability using PIC simulations and establish the con- 
ditions under which it is present. In order to understand 
the saturation mechanisms of the instability, we sepa- 
rate our study in three parts, presented in ^ ^ and 
Q In ^ we show the main non-linear properties of 
the CRCD waves, focusing on their intrinsic saturation 
mechanism in the presence of constant CR current, Jcr- 
We do this using an one-dimensional, analytic model, 
and check our results with one-dimensional PIC simula- 
tions. We calculate a non-linear dispersion relation that 
includes the time evolution of the phase velocity of the 



As mentioned by |Ben| ( [2004|>, the polarization of the waves will 
be right-handed when Bq and Jcr are parallel. In the antiparallel 
case, the polarization is left-handed 



waves. Also, our model quantifies all the plasma motions 
induced by the waves, which is needed to understand the 
wave behavior in the multidimensional context. In ^ we 
present the multidimensional properties of the instability 
using two- and three-dimensional simulations, also with 
constant Jcr- First, we determine the conditions under 
which the CRCD instability can grow without being af- 
fected by plasma filament ation as in the case of ' Niemiec] 
et al.| ([2008). Then, combining the multidimensional sim- 
ulations with our results from ^ we determine the main 
properties of the instability in its non-linear stage. We 
confirm the generation of turbulence as suggested by pre- 
vious MHD studies, and reexamine its main properties. 
We estimate the typical turbulence velocity and length 
scale as a function of the magnetic amplification, finding 
a faster migration to longer wavelengths than predicted 
by MHD simulations. We find that the acceleration of 
background plasma along the direction of motion of the 
CRs causes the intrinsic saLtiiration of the CRCD insta- 
bility at constant Jcr- In q4j we study the effect of the 
back-reaction on the CRs as a second saturation mecha- 
nism for the instability, ^presents our conclusions and 
an application to the case of SNR environments. 

2. CRCD WAVES 
In this section we present the one- dimensional analysis 
of the CRCD waves at constant Jcr, i.e., without con- 
sidering the back-reaction on the CRs. In §2.1 we show 
an analytic, kinetic model for the CRCD waves, valid 
in the non-linear regime. After that, in §2.2, we check 
our model making use of one-dimensional PIC simula- 
tions. We consider a piece of upstream plasma through 
which positively charged CRs ffow, providing the current 
Jcr- We focus on the situation where the initial magnetic 
field Bq^ Jcr^ and the wave vector of the CRCD mode 
are parallel. Local charge neutrality is assumed as ini- 
tial condition, so Ui + Ucr = ^e, where n^, ng, and Ucr 
correspond to the density of ions, electrons, and CRs, 
respectively^. 

2.1. Analytic, kinetic model 
In this section we study the time evolution of CRCD 
waves of different wave vectors /c, which are characterized 
by their growth rate 7 and phase velocity co/k^ where 

k = \k\. The derivation is made for a right-handed po- 
larized wave, and is based on the calculation of the drift 
velocities of plasma components in the presence of a wave 
of arbitrary amplitude. If we know these drift veloc- 
ities and the number densities of the different species, 
we can calculate the total current provided by the back- 
ground plasma as a function of space and time. Adding 
this total plasma current to the constant Jcr contributed 
by the CRs, the time evolution of the wave can be di- 
rectly obtained from the Ampere's and Faraday's laws. 
The details of the calculation are presented in Appendix 
[a) Here we describe its main results, which are summa- 
rized in the dispersion relation given by Equation (A12). 
This dispersion relation assumes a constant 7 and allows 

^ This charge neutrality implies that the background plasma 
must extract a small amount of additional electrons from its sur- 
roundings to compensate the CR charge. 
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oo/k to evolve in time. We will see below that 7 is in- 
deed constant as long as Vd^cr ^ ^a, where va is the 
Alfven velocity of the backgroud plasma. This condi- 
tion not only puts a limit to the validity of the deriva- 
tion, but also sets a saturation criterion for the CRCD 
waves. Also, our derivation is in the low plasma temper- 
ature limit and assumes that va,o ^ {ncr/ni)vd,cr^ where 
va,o is the initial Alfven velocity of the plasma. Using 
two-dimensional PIC simulations, we show below that 
this second condition is actually a requirement for the 
CRCD waves not to be quenched by Weibel-like plasma 
filament at ion. From the real part of Equation (A12), 
we see that the growth rate, 7, is maximized when the 
wavenumber k = kmax = ^TiJcr/B^c^ which corresponds 
to the same wavenum ber of maximum growth found in 
the linear regime (BelTl[2004l i Reville et aL||2QQ6l |Amato 
& Blasi 2008). From the imagmary part, we obtain the 
following differential equation for a; as a function of the 
amplification factor of the waves, / (defined as the ratio 
between the magnitude of the transverse magnetic field, 
Btr, and Bq), 



1 + - 



duo 



^uj[ 1+ 



1 + /^ 



The solution for Equation ([T]) is 

km.n.T.C 



uo = 



0. (1) 



(2) 



which, when va <C c, can be approximated as uo 

This implies that, although in the lin- 



max ^ A 

ear regime the CRCD waves are almost purely growing 
{oj/kmax <^ ^A,Q))i if A g^ts closc to Vd^cn their phase 
velocity can also become comparab le io Vd.cr- 

Taking the real part of Equation (A12) and evaluating 
k = kmax, we obtain that 



at 



^max 



'^max^ 



2 

A,0 



1 



1 



where we have kept terms only to first order in v\/v^ 

and v\l(? . We see that in the regime va ^ "^d.cr, the 
instability grows exponentially with a maximum growth 



rate, 7^ 



that is constant and has the same 



value as obtained in the previous linear studi es (Bell 



2004 
though 



Reville et al 
Equation 



4lf 
IF 



20061 lAmato & Blasi 20081). Even 



shows that our assumption of con- 



stant ^raax is ouly valid when va <C Vd^cr-, it also indi- 
cates that, as Va approaches Vd^cr, the growth rate will be 
substantially reduced, suggesting an intrinsic saturation 
limit for the CRCD waves at va ^ Vd^cr- 

In Appendix [A| we also show that the presence of the 
CRCD waves induces bulk motions of the plasma parti- 
cles both parallel and transverse to Jcr- The parallel mo- 
tion has a velocity Vx^an ^ f'^^A o/^d,cr, while the trans- 
verse motion has a velocity Vtr,an ^ f'^A,o, which always 
points perpendicular to Bfr (and to Jcr)- The parallel 
plasma motion implies that, when va ~ Vd,cr, the entire 
plasma will move at a speed close to Vd,cr, which, from 
the point of view of the plasma, substantially reduces Jcr- 
This reduction in Jcr explains the intrinsic saturation of 
the waves at va ^ yd,cr- The tranverse motion, on the 



other hand, becomes of the order of the Alfven velocity 
of the plasma when the CRCD waves become non-linear 
(/ > 1). We will see below that this increasing trans- 
verse velocity is related to turbulence formation in the 
non-linear regime. 

2.2. One- dimensional Simulations 

We checked our analytical results with one-dimensional 
PIC simulations. We use the PIC code TRISTAN-MP 



( |Buneman|1993l|Spitkovsky|2005D , which can run m one 
two, and three dimensions. In these simulations, like ir 



our analytic model, all plasma properties depend only on 
one spatial direction (x), but both the velocities of the 
particles and the electromagnetic fields keep their three- 
dimensional components. We set up a periodic box that 
contains an initially cold background plasma (with typi- 
cal particle thermal velocity of 10~^c) composed of ions 
and electrons, and a small population of relativistic ions 
(CRs). The driving current is given by the CRs that 
move along x with a mean velocity Vd^cr that we vary be- 
tween runs. These CRs are not allowed to change their 
velocities, as if they had an infinite Lorentz factor, F. 
Such "locked" CRs allow us to study the non-linear evo- 
lution of the instability considering a constant Jcr, i.e., 
eliminating the back-reaction on the CRs. We give elec- 
trons a small velocity along x such that the background 
plasma carries a current —JcrX. This way the net current 
is zero. The initial magnetic field, Bq also points along 
X. Since we want to simulate a situation where ricr <^ ^i, 
having good CR statistics would imply a large number of 
particles per cell. In order to overcome this difficulty, we 
have initialized the same number of macroparticles for 
CRs, ions, and electrons, but modified their charges so 
that Qi = — ge(l — and qcr = — ^e<^, where a = Ucr/rii. 
We change the mass of the particles accordingly in or- 
der to keep the right charge to mass ratios. Particles 
are initially located randomly in the box such that at 
the position of each ion we also have an electron and a 
CR, so the initial charge density is zero. This initial- 
ization is also used in our two- and three-dimensional 
runs. The common numerical parameters for the sim- 
ulations are c/cOp^c ^ 3.3 A (where cOp^c is the electron 
plasma frequency and A is the grid cell size) , ion-electron 
mass ratio m^/me = 100, speed of light c = 0.1125 A/ At 
(where A^ is the time step), and 12.5 particles per species 
per cell. 

2.2.1. Relativistic Regime 

Here we test the non-linear dispersion relation found 
in ^2.1 in the relativistic regime {vd^cr — c), which would 
be appropriate for the upstream medium of a relativis- 
tic shock front. We ran one-dimensional simulations in 
boxes of different sizes L, set up to probe the growth rate 
of different wavelengths, A. As an initial condition we 
used a right-handed, circularly polarized, growing wave 
of amplitude O.l^o, whose fields and particle velocities 
were determined from our analytic model (Appendix IaI). 
We put only one period of the wave in a box, so tnat 
the only other modes that could be excited are shorter 
or equal to L/2. We choose the density of CRs such that 
the corresponding maximum growth rate of the instabil- 
ity, Jmax, is 0.2 cc^ci, SO wc are in the regime where the 
background plasma is well magnetized. Simulations were 
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run for A equal to 0.5, 0.75, 1, 1.25, 1.5, 2, and SXmax- 
We used several values for the initial Alfven velocities, 
'^A,o/'^d,cr, in the range 1/10 to 1/80. This implies that 
the initial gyrotime 27tu;~1 ranges from 189 to 1456 A^, 
so it is resolved with about 20 even when va ^ Vd,cr^ 
and ujp^e/uJc,e ranges from 1 to 8. The results for the 
cases VA,o/'^d,cr = 1/10 and 1/80 are presented in Fig. [l] 
We observe that for L = O.bXmax there is practically no 
growth. For L = 1, 1.25, and l.bXmax we obtain nearly 
the same growth rate, which is very close to the analytic 
Jmax- For longer wavelengths, the growth rate gradu- 
ally decreases. In all our experiments, for L > \maxi 
the exponential growth continues until va becomes close 
to c, which confirms our analytical saturation criterion 
for fixed CRs. At later times we see that, depending on 
L, the amplitude of the wave either oscillates or keeps 
growing but at a much lower rate. 

2.2.2. Non-relativistic Regime 

While individual CRs near a non-relativistic shock 
move at almost the speed of light, on average they move 
with respect to the upstream at a drift velocity, Vd^cr^ 
that is less than c. In order to study this case, we ran 
a series of simulations where, besides not allowing CRs 
to alter their trajectories, we make them drift along x 
at a velocity Vd^cr = 1,0.9,0.8,0.6,0.4, and 0.2c. We do 
this using a box size L = 20Xmax' In this case we do 
not seed the instability with a small amplitude, growing 
wave, as done in §2.2. 1[ Instead, we only put the initial 
magnetic field, Bq^ such that va,o/c = 1/10, forming an 

angle 6 = 5^ with Jcr = JcrX. We use this set-up to show 
that the instability can develop from any kind of noise. 
We tilt Bq by a small angle to inject a small amount of 
magnetic and kinetic energy perpendicular to Jcr, so it 
acts as an initial seed that does not favor any particular 
A (experiments with 6> = were also run, showing no 
difference besides requiring a longer initial time for the 
wave to appear). Other numerical parameters are the 
same as in the relativistic experiments. Fig. [2] shows the 
magnetic energy evolution for the six Vd^cr tested. We 
observe that the growth rate is 7 ^i:^ ^max'^d.cr/ where 
Imax is the maximum growth rate for Vd^cr = c. The 
amplitude at which the exponential growth stops is such 
that Va ^ Vd,cr^ confirming our results from ^2.1 Fig. [S] 
shows the different components of the magnetic field as 
a function of position, x, at different times for Vd^cr = c, 
and Vd^cr = 0.6c. We can see that in both cases the insta- 
bility appears as a right-handed polarized wave and at 
an initial wavelength A ^ Xmaxc/^d^cr^ where Xmax is the 
wavelength of maximum growth for Vd^cr = c. After the 
wave reaches saturation, there is a migration into longer 
wavelengths. This migration appears because the modes 
with wavelengths greater than Xmax grow more slowly, 
but stillgrow and saturate at va ^ Vd,cr^ as can be seen 
in Fig. IT] for the relativistic regime {vd^cr = c). It means 
that, as the instability reaches va ^ Vd,cri the spectrum 
of the waves gradually receives more contribution from 
wavelengths longer than Xmax- 

2.2.3. Background plasma motion 

As we saw in §2.1[ CRCD waves induce plasma motions 
both parallel and perpendicular to x. The left panel in 
Fig. [4| shows the mean velocity of plasma particles along 



X (i.e., parallel to Jcr), and the analytic estimate for 
this velocity, Vx^an ^ f'^'^A o/^d^cr- This velocity comes 
from the E x B drift of background particles, which in 
a well magnetized plasma {jmax <^ ^c,i) is much larger 
than other plasma drifts (see Appendix [A|). The veloc- 
ities in Fig. |4] are computed for two simulations from 
the right panel of Fig. [ij Vd,cr = c with A = 1.5Amaa; 
and Vd^cr = c/2 with A = oXmax' These wavelengths cor- 
respond to the fastest growing waves for the two Vd^cr- 
The right panel in Fig. |4j shows the mean magnitude of 
the transverse velocity of plasma particles for the same 
simulations, and the analytic estimate, Vtr,an ^ /'^a,o- 
Considering that these two cases keep growing exponen- 
tially until tjmax ^ 7 and 12, respectively (see Fig. [T]), 
we see that, in the va <C Vd^cr regime, our analytical es- 
timates for both longitudinal and transverse motions are 
in good agreement with our numerical results. 

3. MULTIDIMENSIONAL EFFECTS 

So far we have studied the properties of the CRCD 
waves assuming an ideal one-dimensional geometry and 
constant CR current. In this section, we relax the first of 
these conditions and use two- and three-dimensional PIC 
simulations to study the CRCD instability, still keeping 
Jcr constant. We identify two main differences with re- 
spect to the one-dimensional case. The first has to do 
with the possibility of plasma filamentation that hap- 
pens before the CRCD instability sets in, as suggested 
by previous works (Nie miec et al.||2Q 08). We will see be- 
low that this filamentation does not occur if the plasma 
is sufficiently magneitzed, Vd^cri'^cr/^i) <^ '^a,o- The sec- 
ond multidimensional effect is the interference between 
CRCD waves generated in different regions of space. 
Since typically the instability starts from random noise, 
different regions will give rise to CRCD waves which in 
general are out of phase with each other. During the 
non-linear stage, this non-coherence makes the transverse 
plasma motions from adjacent regions interfere with each 
other, giving rise to density fluctuations and turbulence 
in the plasma. 

3.1. Magnetization requirement 



Motivated by previous PIC studies by [Niemiec et al. 
(2008), we studied the possibility of an initial plasma 
filamentation that could suppress the formation of the 
CRCD instability. We ran a series of high space resolu- 
tion (c/cjp^e = 10 A) two-dimensional simulations whose 
numerical parameters and results are described in Table 
[1] All our two-dimensional simulations are set up in the 
x — y plane, with Jcr and Bq parallel to the x axis. Also, 
as in some of our one-dimensional simulations, there is 
a small component of Bq pointing along i, working as a 
seed for the instability. 

We identified three regimes, represented in Figs. [5j 
[6] and [71 Fig. [5] shows the plasma density and three 
components of the magnetic field for a simulation with 
'^d,cr(^cr/^i)/'^A,o = 4 (ruu M5 in Table IT]) at two times 
timax = 3 and 11. Even though some (JRCD field is 
observed, especially in By^ the dominant instability cor- 
responds to a transverse filamentation that appears ini- 
tially on the scale of ~ 10 times the electron skin depth. 
As time goes on, the filaments merge, creating promi- 
nent holes in the plasma that preclude the growth of 
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Fig. 1. — Magnetic energy growth, B^^/Bq, is plotted as a function of time, t, for one-dimensional runs with constant Jcr- All the runs 
(except for the one represented by the blue, dot-dashed line in the right panel) are relativist ic (vd^cr — c). In each case we use as a seed 
a single CRCD wave whose wavelength. A, is the size of the simulation box, and the initial amplitude is O.IBq. In all the relativistic cases 
the density of CRs is such that the theoretical maximum growth rate, ^max = 0.2uJc,i, so 7max <C oOci- ^Ao/c = l/lO for the runs in 
the left panel and 1/80 for the ones in the right panel. The green and blue dotted lines correspond to modes with wavelength A =0.5 and 
0.75Amax, the black solid line to A = Xmax-, and the black, blue, green, and red dashed lines correspond to A = 1.25, 1.5,2, and 3Amax, 
where Xmax is the wavelength of the theoretically determined fastest growing mode. The black dotted line represents B'^ in all cases. The 
blue dot-dashed line in the right panel corresponds to a non-relativistic counterpart of the va,o/c = 1/80, A = 1.5Amacc simulation (which 
is the one that grows the fastest), i.e., they have the same parameters except Vd,cr/c = 0.5 and A = SXmax- The results are consistent with 
the theoretical 7max and Xmax, and with the intrinsic saturation criterion, va ~ '^d,cr- 
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Fig. 2. B'i^/Bl is plotted as a function of time, for exper- 
iments similar to the ones depicted in Fig. Jl] but for va,q/c = 
1/10 and a box of 20Amaa;, where Xmax is tne wavelength of the 
theoretically determined fastest growing mode. The instability is 
seeded by tilting Bq by ~ 5° with respect to Jcr- A series of Vd^cr 
is tested: Vd^cr/c = 1 (solid black), 0.9 (solid green), 0.8 (solid red), 
0.6 (dotted black), 0.4 (dotted green), and 0.2 (dotted red), ^max 
is the maximum theoretical growth rate for the case Vd^cr/c- = 1. 
The results are consistent with the theoretical jmax and Xmax-, 
and with the intrinsic saturation criterion, va ~ Vd,cr- 

the instability. Fig. [7| on the other hand, shows the 
same quantities for a simulation where the relative num- 
ber of CRs was decreased by a factor of 10 (run M7), 
implying that '?^d,cr(^cr/^i)/'^A,o = 0-4. In this case, a 
CRCD wave of the size of the box does form (we have 
chosen the x-size of the box to be ~ Xmax)- Finally, 
Fig. [g] shows the case '^(i,cr(^cr/^i)/'^A,o = 2 (run M6), 
in which both the CRCD instability and the initial fila- 
mentation coexist (we call these cases "transitional" and 
indicate them with the letter "T" in Table fl]). These 
three examples indicate that the CRCD instability will 
develop as long as Vd,cr{'^cr/'f^i) ^ '^a,o, which is equiv- 
alent to having a well magnetized plasma in the sense 
that 7max <^ ^c,i (see Appendix [A|). As shown in Table 
[T] we tested the dependence of this criterion on both the 
magnetization of the plasma (using ujp^e/^c,e = 3.15 and 



Vd.cr = c (t7„,,, = 9) Vd_„ = -Sc (ty^ax = 14) 




Fig. 3. — Transverse components of the magnetic field. By (red 
line) and (green line), are plotted as a function of distance, 
X, at two different times for two of the runs described in Fig. [3] 
i'Vd,cr/c = 1 and 0.6). The black line represents Bx- The two 
left plots show the case Vd,cr/c = 1 at tjmax = 9 and 50. The 
two right plots show the case Vd,cr/c = 0.6 at tjmax = 14 and 50. 
The results are consistent with the instability appearing initially as 
right-handed circularly polarized wave with a preferred wavelength 
~ Xmax(vd,cr/c), and with a migration into longer wavelengths as 
the instability grows. 

31.5) and the mass ratio, m^/me (using m^/me = 10 and 
100). We see no difference in our results except that the 
runs with m^/mg = 100 require a slightly higher value of 
the ratio Vd^cri'^cr /^i) / va,o (a factor of 2 larger) for the 
transverse filamentation to dominate, but the qualitative 
criterion remains the same. Also, varying mi/rrie allows 
us to determine the physical length scale, Xfn^ at which 
the transverse filaments appear. This scale shows no de- 
pendence on mi/rrie and corresponds to ~ 10 times the 
electron skin depth. Another mi/rrie = 100 simulation 
was run with zero initial magnetic field (M8), showing 
that the filaments appear at practically the same scale 
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Fig. 4. — The mean velocity along Vx (left panel), and the mean magnitude of the transverse velocity, vtr (right panel), as a function 
of time. Two simulations described in the right panel of Fig. ^(-^^,0/^ = 1/80) are displayed, namely, a simulation with A = l.SAmax for 
'^d,cr/c = 1 (solid lines) and a simulation with A = SXmax for Vd,cr/c = 0.5 (dashed lines). The velocities are represented using red lines. 



For comparison, our analytic estimates for the plasma velocities, Vx 



f '^A Q/'^d,cr and Vtr, an 



seen in the rigth panel of Fig. ^ saturation for the Vd^cr/c = 1 and 0.5 runs happens at tjn 
estimates for the plasma velocities are, therefore, consistent with our numerical results. 



- f'^A,o^ are shown in black lines. As 
7 and 12, respectively. Our analytic 



as in the finite magnetic field case, suggesting a similarity 
between this filamentation and the Weibel instability. Fi- 
nally, these results were also tested for a non-relativistic 
case, Vd^cr = c/2, obtaining the same conclusions. 

Thus, the CRCD instability will grow if the condition 
Vd,cr{ncr/ni) < VA,o (or, equivalently, jrnax < ^c,i) is 
satisfied. For co mparison, the smallest Vd,cr{'^cr/ni)/vA,o 
factor used by iNiemiec et al.| (|2008|) is 1.31, which is 
close to the regime where the transverse filaments ap- 
pear. This fact would explain their plasma filamenta- 
tion, which may have suppressed the appearance of the 
CRCD instability. 

3.2. Multidimensional Evolution 

We studied CRCD instability with a series of two- and 
three-dimensional simulations whose numerical parame- 
ters are summarized in Table [2] As we will see below, 
when multidimensional effects are considered, the dom- 
inant wavelength of the instability, A^^, is initially equal 
to Xmax but then rapidly grows as the field is amplified. 
This can make A^^ equal to the size of the box L be- 
fore the instability reaches saturation, which can make 
sufficiently large three-dimensional simulations challeng- 
ing. We discuss here the results of our three-dimensional 
runs, and check them in Appendix [B] with large two- 
dimensional simulations, for which A^^ is always signifi- 
cantly smaller than the size of the box. 

3.2.1. Three-dimensional Simulations 

In this section we present the results of three three- 
dimensional simulations that test saturation in the non- 
relativistic and relativistic regimes, and the dependence 
of the amplification on the initial magnetic field and the 
CR drift velocity. Two of the simulations have the same 
va,o/c = 1/40, but Vd^cr = c and c/2 (runs II and 12 in 
Table |2|. The third simulation has va,o/c =1/10 and 
^d,cr = c (run 13 in Table |2|. As in all our simulations so 
far, Jcr and Bq point along x (apart from a small compo- 
nent of the magnetic field along z of magnitude Bq/IO)^ 
and the back-reaction on the CRs is not included. The 
rest of the numerical parameters are specified in Table [2] 



The evolution of the plasma density and the three com- 
ponents of the magnetic field for simulation II can be 
seen m Figs. H [lOl and [11] These fi gures show two 
slices of the simulation box. One is longitudinal and cor- 
responds to the plane z/X^ax = 8 (top panels), and the 
other is transverse and corresponds to x/Xmax = 8 (bot- 
tom panels), where Xmax = L/16. The magnetic energy 
evolution for the same run is depicted in Fig. [Ts] 

Fig. [s] shows the early moments of the instability {t = 
47"^^). The longitudinal slice shows how the CRCD 
waves form independently in different regions of the box. 
This is also seen in the transverse slice, which shows how 
the phases of the waves differ between different points of 
the plane x/Xmax = 8- 

Fig. [9] shows the beginning of the non-linear regime 
(/ ~ 1), which corresponds to tjmax = 7. In this case, 
the phases of the waves are transversely more correlated 
compared to Fig. [Sj as can be seen in the plots of the 
transverse slice for By and B^. This increased spatial 
correlation indicates that, until this moment, the adja- 
cent waves were merging without significantly interfer- 
ing with each other. Also, Fig. |9] shows the appear- 
ance of prominent density fiuctuations (An^/n^^ 10, 
where Up is the plasma density). As we saw in q2j when 
the CRCD waves are in the exponential growth regime 
{va <C Vd^cr)^ the background plasma will move trans- 
versely at Vtr,an ^ '^A,o/- So, whcu the instability gets 
non-linear, the transverse velocity of the plasma becomes 
close to vA' In a low temperature regime, this veloc- 
ity corresponds to the magnetosonic sound speed of the 
plasma. Thus, as soon as Bfr ^ Bq^ the transverse mo- 
tions will produce moderate shocks, giving rise to signif- 
icant density fiuctuations in the plasma. At this point 
the transverse plasma motions develop into isotropic tur- 
bulence with velocities of the order of va,o/- What hap- 
pens after the density fiuctuations appear can be seen 



in Fig. [10 {t = 9^max)' "^he longitudinal slice shows 
how the magnetic fiuctuations get distorte d and increase 
rapidly in size. As already mentioned in ^ 2.2.2[ even in 
one-dimensional geometry the instability is expected to 
evolve into wavelengths longer than Xmax- However, in 
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Fig. 5. — Plasma density, rip, and the three components of the magnetic field: B^-, By, and Bz for the two-dimensional simulation M5 
(see Table Q at two different times: tjrnax = 3 (top panels) and 11 (bottom panels), for a box size of lOXmax x 20Amaa;, where ^max and 
Amax are growth rate and wavelength of the fastest growing mode. The arrows in the Bx panels represent the direction of the magnetic 
field projected on the x — y plane. In this simulation Vd^cripcr /Tii) = 4:Va,o- The Weibel-like plasma filamentation on scales of ~ 10 electron 
skin depths can be seen in the plasma density and Bz plots of the top panels. The bottom panels show how the density fiuctuations grow 
both in size and amplitude, suppressing the growth of the CRCD waves. 



plasma density, n, 




Fig. 6. — Same as in Fig. [s] but for simulation M6 of Table ^ The only difference with run M5 is that the density of CRs is reduced 
by half, so cri'^cr/Tii) = '2va,o (which implies an mcrease m A^ax 

by a factor of 2). In this case the CRCD waves and the Weibel-like 
filamentation coexist. We call this situation "transitional" and represent it by the letter "T" in Table \\\ 



a multidimensional set-up, the evolution into magnetic 
fluctuations of larger size gets accelerated after the ap- 
pearance of the density fluctuations and turbule nce in 
the plasma. We will quantify this migration in §3.2.2[ 
The transverse slice of Fig. [lO] shows how the under- 
dense regions (or holes) have merged and increased their 
size with respect to Fig. [9) It is also interesting to see 
from Fig. [12] how the holes are separated by "plasma 
walls" through which the transverse magnetic field re- 
verses direction. We see that the magnetic field has a 



clockwise orientation around the holes, which is consis- 
tent with the presence of right-handed waves producing 
the expansion of the holes. 



Finally, Fig. 



11 



{t = ll7~Lj shows essentially no 
difference with respect to Fig. |1Q| besides the growth of 
the size of both the magnetic fluctuations and the plasma 
holes, which at this point are close to L. 

Fig. 13 shows the magnetic energy evolution for the 
three-dimensional simulations. We can see that, in the 
three cases, the departure from the exponential growth 
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TABLE 1 

Numerical parameters and results for the "M" simulations. 



Run 


rrii/me 


Ucr/rii 


VA,0/Vd,cr 






fmax 

(At) 


^c,i/ ^max 


T/ ^max 


Lx X Ly 


^max 

(A) 


A/ Xmax 




Ml 


10 


0.4 


1/10 


C 


3.15 


352 


0.5 




1024x2048 


99 




8.5 


M2 


10 


0.2 


1/10 


C 


3.15 


703 


1 


T 


1024x2048 


199 


T 


T 


M3 


10 


0.04 


1/10 


c 


3.15 


3,516 


5 


0.82 


1024x2048 


994 


1.03 




M4 


10 


0.2 


1/100 


c 


31.5 


703 


0.1 




1024x2048 


20 




9.8 


M5 


10 


0.04 


1/100 


c 


31.5 


3,516 


0.5 




1024x2048 


99 




8.5 


M6 


10 


0.02 


1/100 


c 


31.5 


7,031 


1 


T 


1024x2048 


199 


T 


T 


M7 


10 


0.004 


1/100 


c 


31.5 


35,156 


5 


0.43 


1024x2048 


994 


1.03 




M8 


100 


0.1264 





c 


oo 









3500x1024 






13 


M9 


100 


0.1264 


1/63.24 


c 


6.31 


3,516 


0.25 




3500x1024 


157 




9.3 


MIO 


100 


0.0632 


1/63.24 


c 


6.31 


7,031 


0.5 


T 


3500x1024 


314 


T 


T 


Mil 


100 


0.0316 


1/63.24 


c 


6.31 


14,063 


1 


0.43 


3500x1024 


629 


0.7 




M12 


10 


0.253 


1/15.81 


c/2 


10 


1,111 


0.5 




2048x2048 


95 




8.5 


M13 


10 


0.1265 


1/15.81 


c/2 


10 


2,222 


1 


T 


2048x2048 


190 


T 


T 


MM 


10 


0.0632 


1/15.81 


c/2 


10 


4,444 


2 


0.69 


2048x2048 


379 


0.8 





Note. — We list the ion-electron mass ratio m^/me, the CR-ion number density ratio ricr/ni, the ratio between the initial Alfven 
velocity of the plasma and the drift velocity of the CRs VA,o/vd,cr, the CR drift velocity Vd^cr-, the ratio between the plasma and the 
cyclotron frequencies of electrons a;p,e/^c,e, the inverse of the maximum theoretical growth rate 7macc in units of the time step At, the 
ratio between the cyclotron frequency of ions and "ymax-, the measured growth rate 7 divided by ^max-, the dimensions of the simulation box 
Lx X Ly in units of of the grid cell size A, the theoretical wavelength of maximum growth Xmax in units of A, the measured wavelength 
of maximum growth A divided by Xmax-, and the size of the transverse Weibel-like filaments Xfu in units of the electron skin depth As,e 
(= c/uop^e)- The symbol "-" and the letter "T" indicate the cases with no CRCD field growth and the ones with a combination of CRCD 
growth and Weibel-like filamentation, respectively. In all the runs the number of particles per cell, Nppc, is 4. 



occurs shortly after the wave becomes non-hnear (which 
coincides with the generation of significant density fiuc- 
tuations and turbulent motions in the plasma). We also 
see that the final saturation satisfies the va ~ Vd^cr condi- 
tion, which suggests that, when multidimensional effects 
are included, the intrinsic saturation of the CRCD insta- 
bility is still given by the va ~ Vd^cr criterion. Unfortu- 
nately, in our three-dimensional simulations, this satura- 
tion happens when the dominant sizes of the holes and 
magnetic fiuctuations have already become close to ~ L. 
Saturation still happens at va ~ ^d,cr because, when 



Xd = the three-dimensional simulations behave more 
like the one-dimensional simulations presented in §2.2.1 



in the sense that there is only one dominant mode that 
saturates at va ~ '^d,cr- After Xd = the density fluctu- 
ations almost disappear and the turbulent motions trans- 
form into more coherent transverse plasma motions. In 
any case, the va ^ Vd^cr saturation criterion is confirmed 
by two-dimensional simulations presented in Appendix [B] 
for which Xd is always smaller than L. 

Fig. 13 also shows that, in the three-dimensional simu- 
lations, the magnitude of the magnetic component along 
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TABLE 2 

Numerical parameters of the "I" simulations. 



Run 


X Ly{xL^)/A^ 


Xs,e/A 


^max / A 


rrii/me 




'^d,cr 






^c,i/ ^max 


Nppc 


11 


5123 


0.5 


32 


10 


1/40 


C 


12.65 


2083 


3.22 


4 


12 


5123 


0.5 


32 


10 


1/20 


c/2 


12.65 


2083 


3.22 


4 


13 


512^ 


0.5 


32 


10 


1/10 


c 


3.16 


521 


3.22 


4 


14 


2048x4096 


1.25 


124 


10 


1/31.62 


c 


10 


2774 


5 


4 


15 


1024x2048 


2.5 


149 


10 


1/10 


c 


3.17 


949 


3 


4 


16 


1024x2048 


2.5 


149 


10 


1/50 


c 


15.8 


4743 


3 


4 



Note. — We list the dimensions of the simulation box Lx x Ly(xLz) , the electron skin depth As,e (= c/up^e), and the theoretical 
wavelength of maximum growth Xmax in units of A. Also, we list the ion-electron mass ratio rrii/me, the ratio between the initial Alfven 
velocity of the plasma and the drift velocity of the CRs ^'A,o/'^d,cr5 the CRs drift velocity Vd^cr-, the ratio between the plasma and the 
cyclotron frequencies of electrons ujp,e/^c,e-, the inverse of the maximum theoretical growth rate ^max in units of the time step At, the 
ratio between the cyclotron frequency of ions and jmax-, and the number of particle per cell Nppc. 



plasma density, n 




Fig. 8. — Plasma density, np, normalized using the mean plasma density, 

'^p,mean: and the th ree components of the magnetic field! Bx^ 
By, and 5^, normalized in terms of ^o, for the three-dimensional simulation II described in Table[2] The top and bottom panels correspond 
to longitudinal (z/Xmax = 8) and transverse {x/Xmax = 8) slices of the simulation box, respectively. At this moment (tjrnax = 4), the 
CRCD waves are still in the linear regime, growing exponentially at a rate ~ jmax (see Fig. \13\ and at a preferred wavelength ~ Xmax- 
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Fig. 9. — Same as in Fig. [8|but for tjmax = 7 (when Btr ~ ^o)- Prominent holes are starting to form in the plasma {Aup/up ~ 10) on 
scales of ~ Xmax- Also the URCD waves are starting to get distorted by the turbulent motions in the plasma. Also, the transverse slices 
show how the transverse spacial correlation of the waves have increased compared to the t^max = 4 case (depicted in Fig. [sj. 

X is comparable to the transverse one, suggesting a rather 
isotropic orientation of the CRCD field. 



3.2.2. Migration into longer wavelengths 

Even though migration to longer wavelengths is al- 
ready observed in one-dimensional simulations, it be- 
comes faster when multidimensional effects are consid- 
ered. In this section we propose a semi-analytic model 
that quantifies this migration in terms of t he am plifi- 
cation factor of the field, /. As we saw in the 
motions associated with the turbulence tend to distort 



the CRCD waves, producing a damping of the shortest 
wavelength modes. Thus, the dominant wavelength, A^^, 
will correspond to the fastest growing mode that can be 
amplified without being strongly affected by the turbu- 
lence. Considering that a CRCD wave of wavelength 
A grows in a time scale comparable t o the inverse of 



its growth rate 7"^ (A) from Equation (A12), and that 
the turbulence will kill it in a time scale comparable to 
X/vturb^ where Vturb is the typical turbulent velocity, then 
Xd will be such that j~^{Xd) ^ X/vfurb- Since the turbu- 
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Fig. 10. — Same as in Fig. ^but for tj max — 9 (Btr ^ ^0, see black line in Fig. |13| ). The qualitative features of the turbulence are the 
same than in Fig. [o] except for an increase in the size of the plasma holes and dommant CRCD wavelength. 




plasma density, n 




Fig. 12. — The density plot presented in Fig. [Toj but with 
overplotted arrows showing the magnetic field projection on the 
z — y plane. The clock-wise orientation of the magnetic field lines 
around the plasma holes shows the presence of CRCD waves driving 
the turbulence. 



lence is due to the transverse plasma motions produced 
by non-coherent CRCD waves, then Vturb must be com- 
parable to the transverse velocity of the waves, which we 
already determined to be /va^o- So, A^^ will be such that 
li^d) ^ /3vA,of / where /3 is an unknown constant 
that quantifies the relative importance of the two time 




Fig. 13. — The transverse (solid) and longitudinal (dotted) com- 
ponents of the magnetic energy as a function of time for three- 
dimensional simulations 15 (black), 16 (red), and 17 (green) of Ta- 
ble [2] for which (vA,o/vd,cr,Vd,cr/c) = (1/40,1), (1/20,0.5), and 
(1/10,1), respectively. Time is normalized in terms of the jmax 
of each simulation. In all the runs, the departure from expo- 
nential growth occurs after Btr ~ ^0, but saturation happens at 

VA ~ Vd,cr- 

scales. If we get ^{Xd) from the real part of Equation 



(A12), we can obtain P by fitting the evolution of A^^ 
in our three-dimensional simulations, obtaining (3^2. 
This way we find A^^ as a function of /, 



max [(//3)' + l]/2, 



(4) 



which is intended to be valid after the turbulence be- 



comes significant (/ > 3). In Fig. 14 we show a com- 
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Fig. 14. — Evolution of the dominant wavelength, A^f , as a func- 
tion of the amplification factor, /, for three-dimensional simula- 
tions II (dot-dashed), 12 (dashed), and 13 (dotted). For compar- 
ison, our semi-analytical formula, = Xmax((f f^)"^ + l)/2, is 
shown as solid line. 

parison between this formula and the evolution of A^^ as 
a function of / for our three-dimensional simulations. 
We computed A^^ by performing Fourier transforms of 
By and Bz along lines of constant y and z coordinates, 
and then finding the mean wavelength of the peak of 
the Fourier transform. As the dominant wavelengths ap- 
proach L (= WXmax)^ the determination of A^^ becomes 
quite noisy. Due to this reason, we have plotted our sim- 
ulation results only until A^^ = L/2. We see from Fig. 
14 that Equation (|4| appears to provide an acceptable 
nt for the evolution of A^^. Our result shows a growth 
of Xd substantially faster than t he di r ect pr oportionality 
between A^^ and / suggested by Bell (2004). We will see 
below that, when the back-reaction on the CRs is con- 
sidered, this difference has important implications to the 
saturation of the CRCD instability. 

4. BACK-REACTION ON COSMIC RAYS 

In this section we use one- and three-dimensional sim- 
ulations to study the effect of the dynamic evolution of 
the CRs on the saturation of the CRCD instability. We 
will concentrate on the case of a beam of monoenergetic 
CRs drifting at half the speed of light {vd^cr = 0.5c) in 
the X direction, which is parallel to Bq (except for a small 
magnetic component along z). This drift is obtained by 
sampling the CR velocities from an isotropic, monoener- 
getic momentum distribution, but only keeping the ve- 
locities in the positive x direction. This choice for the 
CR momentum distribution has a direct application to 
the most energetic CRs that propagate in the upstream 
medium of SNR shocks (see 

The red and green curves in Fig. 15 represent the 
magnetic energy evolution for two one-dimensional sim- 
ulations that only differ in their Lorentz factors, F, taken 
to be 20 and 40 (the rest of the numerical parame- 
ters are specified in the caption of Fig. fT5|. These 
simulations saturate at / ~ 7.9 and 12.6, when the 
Larmor radius of CRs is close to the dominant wave- 
length of the instability {Rh^cr/Xd = 1.05 and 1.12 for 
two runs). To obtain this result we use that initially 
RL,cr/Xmax = {ricr / ni){vd,cr / c){c/va,o)'^T / Att , and con- 
sider that at saturation the dominant wavelength has 
grown {Xd/Xmax ~ 2.9 and 3.4, respectively) and CRs 
have lost part of their energy (the mean F of CRs is 18.9 
and 37.6, respectively). These results are confirmed by 



and normalized in terms of B^, are plot- 



FlG. 15. 

ted as a function of time (normalized using jmax) for one- and 
three-dimensional runs in order to study the effect on the mag- 
netic energy evolution due to the back-reaction on the CR. In all 
cases CRs are monoenergetic and have a semi-isotropic momentum 
distribution such that Vd,cr = c/2. The two one-dimensional simu- 
lations have numerical parameters: ricrl'^i — 0.005, nr 

/c = 0.5, 

'^A,{)hd,cr = 1/40, Xmax = 2048 A, L/Xmax = 15, mi/uie = 100, 

^p,e/^c,e = 8, = 231,786 At, uJci/jmax = 10, and 

c/ujp^e = 3.3. Their CR Lorentz factor F is 20 (solid, red line) 
and 40 (solid, green line), respectively. The three-dimensional sim- 
ulation, whose transverse magnetic energy is represented by the 
solid, blue line, has the same parameters as run 12 in Table [2] but 
with r = 30. The dotted, blue line represents the longitudinaTmag- 
netic energy (B'^/Bq). The dashed, black line shows the constant 
magnetic energy along x for the two one-dimensional simulations. 

a three-dimensional simulation whose magnetic energy 
evolution is represented by the blue lines in Fig. TE\ The 



numerical parameters of this three-dimensional simula- 
tion are the same as in run 12 (see Table [2]) but includes 
the back-reaction on the CRs, whose F = 30. If we con- 
sider that at saturation Xd/Xmax = 4 and the mean F of 
the CRs is 27.6, we obtain that the CR deflection satu- 
rates the instability when RL,cr/^d = 1-35. Note that in 
this simulation A^^ is always a factor of 4 smaller than the 
size of the box, so the saturation is not affected by box 
effects. For the semi-isotropic distribution of monoener- 
getic CRs presented here, we find that the saturation due 
to CR back-reaction will happen when Rl^cv ^ Xd- Al- 
though this result is valid for our particular choice of CR 
momentum distribution, we expect that in general the 
saturation of the CRCD instability will be determined 
either by the intrinsic limit va ~ ^d,cr, or by the strong 
CR defiection when Rl^cv ^ Xd. As will be discussed in 
31 achieving va ^ Vd,cr requires a very high CR energy 
ensity, a condition that is not expected for non-relatistic 
shocks environments. 

Also, our simulations show that, at saturation, many 
CRs have negative x velocity. This suggests that, besides 
the field amplification, the CRCD instability can provide 
an efficient scattering mechanism for CRs upstream of 
shocks. 

5. DISCUSSION AND CONCLUSIONS 

Using fully kinetic PIC simulations, we confirmed the 
existence of the CRCD instability predicted by Bell 
(2004). Combining one-, two-, and three-dimensional 
simulations with an analytic, kinetic model we studied 
the non-linear properties of the instability and its possi- 
ble saturation mechanisms. 

In the first part, we studied non-linear CRCD waves 
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under idealized conditions, namely: i) ignoring multidi- 
mensional effects, and ii) assuming a constant CR cur- 
rent without back-reaction on the CRs. We confirm that 
the CRCD waves can grow exponentially at the wave- 
lengths and rates predicted by the an alytic dispersion 



relation (Beh 2004; Reville et al. 2006 Amato & Blasi 



[2008 ) . We find that the exponential growth can continue 
mto the very nonlinear regime, until the Alfven velocity 
in the amplified field is comparable to the CR drift veloc- 
ity, va ~ Vd^cr- This saturation is due to plasma accel- 
eration along the direction of motion of the CRs, which 
reduces the CR current observed by the plasma particles. 
The plasma moves at the velocity ^ f^^A o/'^d,cr^ where 
/ is the amplification factor of the field (/ = Btr/Bo). At 
saturation, when va ^ '^A,o/ ^ "^d.cn the plasma moves 
together with CRs, decreasing the net driving current. 
The waves also induce transverse plasma motions with 
velocities ~ va,o/- These motions generate plasma tur- 
bulence when multidimensional effects are included. 

In the second part, we considered more realistic con- 
ditions by including the multidimensional effects using 
two- and three-dimensional simulations with constant 
Jcr- Our main results are: 

i)ln the linear regime, if the plasma is well magnetized 
{imax < ^c,i, or, equivalently, Vd,cr{ncr/ni) <C va,o), the 
CRCD waves grow at the rate and preferred wavelength 
close to the ones obtained in the one-dimensional anal- 
ysis. If this condition is not met, Weibel-like filaments 
form in the plasma, supressing the appearance of the 
waves. In this case, the streaming CRs can still amplify 
the magnetic field to non-linear values, but at a rate sig- 



nificantly lo wer than that of the CRCD waves (Niemiec 



et al. 2008). This regime, however, might be relevant 



to the upstream medium o f relativistic shocks in GRBs, 
where ricr could exceed (Couch et al.| |2008 ). 

ii) In the non-linear CRCD regime, the transverse 
plasma motions associated with the instability create 
significant density fluctuations and turbulence in the 
plasma. These turbulent motions suppress the growth 
of the shortest CRCD waves, producing a fast evolution 
into longer wavelengths that can be approximated by 
Ad ~ Xmaxiif/^)'^ + Also, even though the field 

will continue to be amplified until va ~ Vd,cr^ the non- 
linear growth will be slower than in the linear regime. 

In the third part, we include the back-reaction on the 
CRs and find that the CR deflection by the amplified 
field constitutes another possible saturation mechanism. 
We tested this effect for a semi-isotropic distribution of 
monoenergetic CRs propagating at Vd,cr = 0.5c with 
respect to the upstream medium, which would be ap- 
propriate for the most energetic CRs that escape from 
SNRs. We find that the field is amplified until the Lar- 
mor radii of the CRs becomes approximately equal to 
the size of the dominant magnetic fluctuations. When 
that happens, the CRs get strongly deflected by the mag- 
netic field, which decreases their current and stops the 
growth of the field. Ignoring the migration to longer 
wavelengths, for a generic CR momentum distribution 
with Vd,cr < c, saturation due to CR deflection will hap- 
pen when / ^ {Tcr/4:7r){c^/v\^Q){ncr/ni){vd,cr/c), where 
Per is the typical Lorentz factor of current-carrying CR. 
On the other hand, saturation due to plasma accelera- 
tion to Vd^cr velocity occurs when / ^ Vd,cr/vA,0' Thus, 



the deflection of CRs will dominate if the CR energy den- 
sity is such that Terser < "^T^^iVA^o / If the migration 
to longer wavelengths is included, saturation due to CR 
deflection would happen at even smaller magnetic ampli- 
fication. 

In the upstream medium of SNR forward shocks, we 
expect the CR energy density to be low enough so that 
the maximum CRCD amplification is determined by the 
back-reaction on the CRs. In order to make an estimate 
of typical magnetic amplification in these environments, 
let us consider a piece of upstream whose distance from 
the shock is such that it can only feel the most ener- 
getic CRs that escape from the remnant. We use only 
the most energetic particles because they are the only 
ones whose Larmor radii are much larger than the typi- 
cal wavelength of the CRCD waves, which is an essential 
condition of the instability. Also, our estimate is based 
on the following assumptions. First, all the escaping par- 
ticles have positive charge, which is reasonable consider- 
ing the much shorter cooling time of electrons compared 
to ions, and that ions are presumably more efficiently in- 
jected into the acceleration process in shocks. Second, 
we assume that the escaping CRs have the same energy, 
Eesc^ which is roughly the minimum energy required for 
them to run away from the remnant. Third, there is a 
fixed ratio, r^esc = ^£;,cr/(p'^s^/2), between the flux of 
CR energy emitted by the shock, Fe^ct-, and the flux of 
energy coming from the upstream medium as seen from 
the frame of the shock, p'u^^/2, where Vsh is the shock ve- 
locity and p is the mass density of the upstream plasma. 
Finally, we assume a plane geometry and that all the 
ions are protons. Under these conditions, the time scale 
of growth of the instability, 7"^^:, is 



^max^^^\ iQiSeV 



lo^km/sec 

Vsh 



0.05 \ ( cm- 

''7esc 



(5) 



years. 



and the initial length scale of maximum growth is 



V lokm/sec 



lo^km/sec 

Vsh 



The ratio jmax/^c,i is 

lmax/^c,i^^'^ X 10" 



0.05 

Vesc 



iQi^eV 
cm- 



(6) 



pc. 



Y iokm/sec \ / iQi^eV 

I VA,0 I \ Eesc 



lo^km/sec 



Vesc 

0.05 



(7) 



which confirms that in the case of SNRs the CRCD in- 
stability will not b e aff ected by the Weibel-like filamen- 
tation studied in ^3.1 Considering the migration into 
longer wavelengths given by Equation (|4|, the amplifica- 
tion factor, /, in SNRs will satisfy 



/[(//3)'+l] 



130 



sec 



(- 



Vsh 



Ve 



0.05 J VlO%m/sec 



(8) 

which, for typical parameters would imply / ^ 10. We 
see that, even though run- away CRs can significantly 
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amplify the ambient magnetic field, the upstream am- 
plification alone is not enough to explain the factors of 
~ 100 inferred from observations of forward shocks in 



et al. 



|2007|. 



young SNRs (jVolketaL 2005 



Also, note that 



Ballet] |2006| 



Uchiyama 



and 



the distance swept by the shock in a time 
0.5pc, which is comparable to the typ- 
ical size of a SNR. This means that the advection of the 
upstream fluid into the shock may happen faster than 
the growth of the field, and may put further restrictions 
on the amplification. 

It has been suggested that a further CRCD amplifi- 
cation could be provided by the current of lower energy 
CRs that are confined closer to the sho ck and mo ve dif- 
fusively at drift velocity Vd^cr ^ Vsh (|Bell||2004|). We 



believe, however, that this possibility requires a more de- 
tailed study. Such lower energy CRs can be magnetized 
in the sense that their Larmor radii are smaller than the 
typical size of magnetic fluctuations, A^^, violating the 
conditions for the CRCD instability. Although on large 
scales these CRs will still produce a current of magni- 
tude ~ ericrVsh parallel to the shock normal, on scales 
of the CRCD wavelength the local CR current may get 
significantly affected by the amplified field because of the 
deflection of CRs. Field amplification may then proceed 
in essentially different way. 

The magnetization of CRs could be less of an issue if 
the wavelength of the instability due to low energy CRs 
is shorter than the CR Larmor radius. Indeed, since 
lower energy CRs are more numerous than the most en- 
ergetic ones, their larger current will generate shorter 
CRCD waves (remember that Xmax — Bqc/Jct)- How- 
ever, the CRCD turbulence generated further upstream 
by the highest energy CRs may modify the condition 
Jcr II and may suppress the growth of the small wave- 
length modes closer to the shock. This suggests that 
other non-linear mechanisms, su ch as the cyclotron reso 



nance of CRs with Alfve n waves (iKulsrud & Pearce 1969 
McKenzie & Volk 1982), may still be important compo 



nents in the amplification of the field. The full effect of 



the low-energy CR contribution needs to be investigated 
using a fully kinetic treatment of CRs that includes their 
reacceleration by the shock, the presence of pre-existing 
turbulence, and the eventual contribution of CR elec- 
trons to the cancelation of ion current. 

Although we have applied our results only to the non- 
relativistic case of SNRs {vd^cr < c), the CRCD insta- 
bility may also play an important role in relativist ic 
shocks in jets and Gamma Ray Bursts, where Vd^cr — c. 
Our simulations show that, at constant CR current, the 
evolution of the instability is the same as in the non- 
relativistic case. In particular, the intrinsic saturation 
criterion due to plasma acceleration is valid, implying 
a maximum magnetic fiel such that va ~ c. However, 
if the back-reaction on the CRs is considered, the non- 
linear evolution of the field and the saturation due to 
CR def lection may be dominated by CR beam filamen- 
tation ( [Riq uelme fc Spitkovsky [2008 ). Also, since in the 
upstream of GRB shocks the aensit y of CRs might b e 
close to the density of upstream ions rtCouch et al. 2008L 
the magnetization requirement {jmax ^ ^c,i) may not be 
satisfied in these environments. In this case, a non-linear 
magnetic amplification is still expected, but through an 
instability that is characterized by Weibel-like filamenta- 
tion of the plasma and whose properties may be different 
to th e CRCD instability described here (Ni emiec et al.| 
|2008[ ). Detailed analysis of the relativistic shock case in 
application to GRBs will be presented elsewhere. 

In conclusion, we have shown that the CRCD instabil- 
ity is a viable mechanism for the non-linear amplification 
of magnetic field upstream of both non-relativistic and 
relativistic shocks, and that it can provide an efficient 
scattering mechanism for CRs in these environments. 
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APPENDIX 
A) CRCD WAVES DISPERSION RELATION 

In this appendix we calculate a dispersion relation for the CRCD waves for the case where Bq^ Jcr^ and the wave 
vector of the electromagnetic mode, /c, are all parallel and point along the x axis. We will separate the fields and 
currents into components that are transverse and parallel to x, and will identify them with the subscripts "tr" (standing 
for "transverse") and "x", respectively. Thus, the magnetic field perpendicular to x will be given by 

Btr = Re{Bo{iy + i)e*(^-^(^))+^^}, (Al) 
which corresponds to a right-handed polarized wave, where the phase 0(t) is an unknown function of time, t, the 
growth rate 7 > is a constant, and Bq is the magnitude of the initial background field ^o- Note that the time is 
chosen so that the wave is in the linear regime for t < 0. 

Then, from the Ampere's and Faraday's laws, we get that the electric field and the current perpendicular to x are 
given by 

■ Bo{u; + ij) ^ 



Re 



kc 



and 



Re 



Bo 



^ Airkc 

where co and co are the time and second time derivative of ( 
we need another expression connecting Jtr with Etr and Btr 



(A2) 



(A3) 



respectively. In order to obtain a dispersion relation. 
We find it by making the following assumptions. First, 
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the thermal velocities of the particles in the background plasma will not give rise to any significant drift velocity. 
Second, we will assume that 7, \d(kx ^ (t)(t)) / dt\ and (d^[cj]/dr)/((i(^-i) [cj]/^^^^-^)) are constant and much smaller 
than \qjB/mjc\ = uJcj, where uJcj, Qj^ and mj are the cyclotron frequency, the charge, and the mass of the j species, 
respectively. (We will see at the end of this appendix that, in order to satisfy the last three conditions, we need 
va,o/c ^ {TicR/Tii){vd^cr/c)^ whcrc va,o IS the initial Alfven velocity of the plasma.) Third, the electric and magnetic 
fields are perpendicular, which is a reasonable assumption in the case of a quasineutral plasma. And finally, Vd^cr ^ ^a, 
where va is the Alfven velocity of the plasma and Vd^cr is the drift velocity of the CRs. Considering this, given Eir 
and Btr^ we can find Jtr as follows. 
If a particle "j" experiences electric and magnetic fields, E and 5, its velocity perpendicular to B has two components, 

/dt = {qj /mj){vgj /c) x B and dvdj/dt = {qj/mj){E + Vdj/c x B). In 
^gj represents the classical gyration around 5, while Vdj corresponds to 



and Vdj^ that satisfy the equations dvg 



the case of constant and uniform E and B 

the drift of the particle, which is Vdj = c{E x B)/B'^. 

When the fields change both in time and space, we can still decompose the velocity perpendicular to the field into 
Vgj -\- Vdj [again, satisfying dvgj/dt = {qj /mj){vgj /c) x B and dvdj/dt = {qj /mj){E -\- Vdj / c x B)]. In this case, the 

space and time variations of B can also produce drift velocities due to the Vgj motion. The space variations will give 
rise to a drift due to the curvature of the magnetic field lines (curvature drift). The curvature drift velocity, however, 
is of the order of the thermal speed of the particles times the ratio between their Larmor radii and the curvature 
radius of the lines. The time variations of the field, on the other hand, can also give rise to drift velocities. To first 
order, these velocities will also be proportional to the thermal speed of the particles times the ratio between the rate 
of change of the field (determined by the quantities 7 and \d{kx + (j){t))/dt\) and oOcj- We will neglect these possible 
drift velocities using our first assumption that the thermal velocities of the particles are low enough not to produce 
any important drift velocity. 

On the other hand, in the case of a non-uniform and time-changing fields, the Vdj velocity is given by the series. 



n=0 



where 



yiO) _ iExB)c . Qjv'J 



X B 



dv 



dt 



(A4) 



(A5) 



forn = 1, 2, 00. We see from Equations (Al), (A2), and (A5) that \v^j\ <C \v^j as long as 7, \d{kx-\-(j){t)) /dt\^ and 

(d^M/dr)/(d(^- 



'j]/dt''^ ^ ^ ) are constant and much smaller than \qjB/ rrij c \ 



which is our second assumption. 

Notice that, even if |^^^]| <C I'^^^jl, the currents produced by these two velocities can be comparably important. This 
is because iJ^^j is independent of rrij and g^, thus it has the same value for ions and electrons (so from now we will 
just drop the subscript "j" and will refer to this velocity as simply v^^). Thus, when considering both species, the 
current produced by ^^^^ to the tiny excess of electrons in the background required to compensate 

the CRs charge, so it will be proportional to ricr- On the other hand, since iT^^j is proportional to rrij/qj^ it will be 
much larger for ions than for electrons. So the corresponding will be proportional to the total density of ions in 
the background, n^, which is typically much larger than ricr- Since v^J^^- for n > 2 will also affect mainly the ions, their 



contribution to the current in the plasma will be much smaller than the one of iT^^- provided that v)!^"^- <C v)^'^ ^^ so 



^(n-l) 



we w ill just neglect them. Thus the currents and can be calculated considering Equations (Al), (A2), and 
( A5 ) , finding that 



J 



and 



kvd. 



^d " 



1 + f 



ri + Re 



1 + 



■1' + 



-r 



Re 



-ly 



(A6) 



(A7) 



where we have defined the field amplification factor / = Btr/Bo. 

This way we have calculated all the currents in the plasma that are perpendicular to 5, but we still have to determine 
the ones that are parallel to B. We do that using our third assumption, E ± B, which implies that 



Ex = —{Etr • Btr)/Bx. 



(A8) 
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plasma density, n 




2 4 5 S 10 12 14 



14 of Table [2] The arrows on the Bx panel show the direction of the magnetic field projected on the x — y plane. This stage shows the 
CRCD waves growing at wavelengths ~ Xmax in their linear stage {Btr <C ^o)- 

Using the Ampere's law and the fact that in a one dimensional problem (V x B)^ = 0, we have that 



dt 



J. 



(0) 

d,x 



J. 



(1) 

djX 



where Jy is the x component of the plasma current parallel to 5, Jy 
to X is just 

Bo' 



^J\l^). (A9) 
Given this, the component of Jy perpendicular 

(AlO) 



Then, using Equations (|Alj), (|A2|, (|A8|), (|A9| and ( |A10[ ) we get that 



Ji 



^tr = Rej ^Jcr^l 



to 







Alike 



■7 



ix) 



l + p 



i{kx-(t>{t))+-ft 



i-iy-z)\. (All) 



Now we have the expressions for all the components of the current perpendicular to f , J^r, so we can use Equation 
(A3) to find the dispersion relation, 



47rJc 
Bo 



22 2 / -1 



27a; ) • 



(A12) 



From this derivation we can also obtain an estimate of the plasma velocities due to the CRCD waves. We know that 



the motion of particles perpendicular to B is dominated by 



J^P /Tier (since v^^^ affects ions and electrons in 



the same way and Ue — rii = ricr)-, and that the motion of particles parallel to B is mainly 



at v\\ = —J\\/ne- Thus, by looking at the expressions for J^^^ and J\\ given by Equations (|A6|), (|A10|), and (|All|), we 



find that the dominant plasma motion will be given by ^f^^ and will imply a velocity of ions and electrons that can be 
decomposed into a component along x, Vx^an ^ f'^'^A o/^d,cr (where the subscript "an" stands for ^^analitic^^)^ and a 

trans verse component, Vtr,an ^ /'^A,o, that is always perpe ndicu lar to Btr (and to x). 

In ^2.1 we use the dispersion relation given by Equation ( |A12 ) to calculate the wavenumber, kmax^ and growth rate, 
Jmax 5 of the fastest growing mode, as well as a; as a function of the amplitude of the wave. We will use these results 
here in order to check the consistency of assuming that 7, \d{kx + (j){t))/dt\, and {d'^[(jo]/dt'^)/{d^'^~^'^[(jo]/dt^'^~^'^) are 
constant and much smaller than uJc,i^ which are necessary for neglecting \v^^\ when n > 2. We know from Equation 
([2]) that when va <C Vd^cr^ ^ ^ kmax'^\/'^d,cr- It means that {d^^\uj\/ dt^) / {d^^~^\uj\/ dt^^~^^) ^ 27 (except when 

n = 1 and / < 1, in which case {d^''^[u;]/dt'')/{d^''-^^[u;]/dt^''-^^) < 7). So, in order to neglect {v^J'^l for n > 2, we 
only require 7 and \d{kx + (j){t))/dt\ to be constant and much smaller than uJci- It is possible to show from Equation 
(|3| that the first condition, which is equivalent to 7mox <^ ^c,i^ where jmax is given by Equation ([3|, is satified if 



given 



by e lectro ns mo ving 



'^A,o/c ^ {TicR/ni){vd^cr/ c)- From Equation (2) we see that, if we approximate dx/dt 



approximately equal to k^naxci'^A^o/cY 



^d,x^ 



d{kx + (j){t) ) / dt becomes 



c('^A,o/c), which is also constant and much smaller than jrnax- This way 



we see that our analytical results are valid if the plasma is well magnetized in the sense t hat 
equivalent to va,o ^ {^CR/^i)'^d,cr- Using two-dimensional PIC simulations, we show in ^3.1 
actually a requirement for the CRCD not to be quenched by the Weible-like filamentation. 

B) MULTIDIMENSIONAL EVOLUTION: TWO-DIMENSIONAL SIMULATIONS 



which is 



that this condition is 



We saw in §3.2.1| that, when multidimensional effects are considered, the dominant wavelength of the CRCD insta- 
bility, Xd, grows according to Equation Q. This makes it numerically expensive to run three-dimensional simulations 
that could amplify the field substantially without making A^^ too close to the size of the simulation box, L. In order to 
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plasma density, 




1 



2 4 6 8 10 12 14 



2 4 6 S 10 12 14 




2 4 5 S 10 12 14 



Fig. 18. — Same as in Fig. 16 but at tj max = 11 (when Btr ^ ^o)- The ampl itud e of the density fluctuations is stih prominent 
(Aup/up ~ 10), but their lengtn scale is larger, compared to the t^max = 9 case (Fig. |17} . The magnetic fluctuations have also grown in 
size, and a clear differentiation between By and Bz, peculiar to two-dimensional runs, can be observed. 

overcome this difficulty, in this section we present the results of two-dimensional simulations whose L is always bigger 
than Xd. Despite some artifacts related to the two-dimensional geometry, these simulations help us confirm the main 
results obtained from the three-dimensional analysis presented above. 

Figures 16, 17, and 18 show the results at three different times {t^max = 6, 9, and 11, respectively) for one of the 
simulations~(run 14 of Table |2|, which corresponds to CRs drifting at the speed of light and without considering their 
back-reaction. We see that initially the instability is produced independently in different regions of the simulation box 
(as seen in Fig. 16). In this linear stage of evolution, the waves produced in adjacent regions of space seem to grow 
without interfering with each other. However, when the waves become non-linear, strong density fluctuations appear 
on scales of a few Xmax (as shown in Fig. 17). The beginning of this stage is shown in Fig. 17 It is also apparent from 



Fig. 18 that the magnetic fluctuations get distorted and evolve into larger scales right after the density fluctuations 
and turbulence form. 

The formation of density fluctuations also affects the growth rate of the instability. Fig. [19] presents the magnetic 
energy evolution for the two-dimensional simulations 15 and 16, whose va,o/c = 1/10 and 1/50, respectively. The rest 
of their numerical parameters are specified in Table [2] We see that, as in the three-dimensional case, the exponential 
growth stops shortly after Btr ^ Bq. After that, the CRCD instability grows at a lower rate, reaching saturation 
when va ~ Vd,cr- This result had already been obtained in the three-dimensional case, but in this case we allow the 
instability to evolve into larger scales as the magnetic field grows. 

In two dimensions, the formation of density fluctuations produces a clear differentiation between the y and z com- 
ponents of the field (as can be seen in Figs. 18 and 19). This is because in the low density regions the plasma cannot 
generate the return current necessary to compensate Jcr- Thus, the uncompensated CR current produces a "toroidal" 
magnetic field around the underdense regions that, in the two-dimensional case, manifests itself as an amplification of 
the out of the plane component of the field, Bz- Even though, as seen in Fig. [18) both the "toroidal" and the CRCD 
field coexist, the two-dimensional simulations can still give us information about the point when the CRCD instability 
stops amplifying the field. 
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